# -*- coding: utf-8 -*-
"""
@Time ： 2022/10/8 12:59
@Auth ： GaoShuai
@File ：ProjectionTest.py
@IDE ：PyCharm
"""
if __name__ == '__main__':
    """Example using the ray transform with circular cone beam geometry."""

    import numpy as np
    import odl

    # Reconstruction space: discretized functions on the cube
    # [-20, 20]^3 with 300 samples per dimension.
    reco_space = odl.uniform_discr(
        min_pt=[-20, -20, -20], max_pt=[20, 20, 20], shape=[300, 300, 300],
        dtype='float32')

    # Make a circular cone beam geometry with flat detector
    # Angles: uniformly spaced, n = 360, min = 0, max = 2 * pi
    angle_partition = odl.uniform_partition(0, 2 * np.pi, 360)
    # Detector: uniformly sampled, n = (512, 512), min = (-30, -30), max = (30, 30)
    detector_partition = odl.uniform_partition([-30, -30], [30, 30], [512, 512])
    geometry = odl.tomo.ConeBeamGeometry(
        angle_partition, detector_partition, src_radius=1000, det_radius=100,
        axis=[1, 0, 0])

    # Ray transform (= forward projection).
    ray_trafo = odl.tomo.RayTransform(reco_space, geometry)

    # Create a discrete Shepp-Logan phantom (modified version)
    phantom = odl.phantom.shepp_logan(reco_space, True)

    # Create projection data by calling the ray transform on the phantom
    proj_data = ray_trafo(phantom)

    # Back-projection can be done by simply calling the adjoint operator on the
    # projection data (or any element in the projection space).
    backproj = ray_trafo.adjoint(proj_data)

    # Shows a slice of the phantom, projections, and reconstruction
    phantom.show(coords=[None, None, 0], title='Phantom, Middle Z Slice')
    proj_data.show(coords=[0, None, None], title='Projection at theta = 0')
    proj_data.show(coords=[None, None, 0], title='Sinogram, Middle Slice')
    backproj.show(coords=[None, None, 0], title='Back-projection, Middle Z Slice',
                  force_show=True)